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Test  on  Marmousi  data  of  velocity  analysis  by 

perturbation 

Zhenyue  Liu  and  Norman  Bleistein 


ABSTRACT 

Velocity  model  determination  on  the  Marmousi  data  is  a  major  challenge  for 
current  velocity  analysis  methods  because  of  the  extreme  complexity  of  the  model. 
A  new  velocity  analysis  method  has  been  developed  and  is  tested  on  the  Marmousi 
data.  This  method  characterizes  a  velocity  distribution  by  the  macro  model  and 
uses  an  analytical  formula  to  update  velocity.  Our  formula  estimates  the  update 
in  velocity  by  computing  a  derivative  function  of  imaged  depths  with  respect  to 
velocity  in  a  general  background  medium  context.  This  formula  is  more  accurate 
than  conventional  ones  based  on  hyperbolic  residual  moveout  when  the  medium 
has  strongly  lateral  velocity  variations.  The  test  results  show  our  velocity  analysis 
method  to  be  superior  for  complex  media  such  as  in  the  Marmousi  model. 


INTRODUCTION 

The  Marmousi  model  is  a  2D  model  with  considerable  complexity  of  structures 
and  a  realistic  distribution  of  reflectors.  Prestack  data  from  it  provide  a  stiff  challenge 
to  methods  for  estimating  velocities  from  seismic  data. 

Strong  velocity  variations  in  the  Marmousi  model  suggest  that  the  method  for 
velocity  estimation  should  not  be  based  on  the  assumption  of  hyperbolic  residual 
moveout.  Liu  and  Bleistein  (1994)  proposed  an  approach  to  migration  velocity  anal¬ 
ysis  that  characterizes  a  velocity  distribution  by  using  the  macro  model — consisting 
of  velocities  and  velocity  interfaces.  The  model  is  determined  by  layer-stripping,  in  a 
top-down  procedure;  the  velocity  of  each  layer  (assumed  constant  or  constant  gradi¬ 
ent)  is  updated  iteratively  by  using  an  analytical  formula;  and  the  velocity  interface 
is  imaged  by  using  the  corrected  velocity.  In  order  to  handle  nonhyperbolic  residual 
moveout,  the  velocity  estimation  in  the  proposed  formula  is  implemented  by  com¬ 
puting  a  derivative  function  of  imaged  depths  with  respect  to  velocity  in  a  general 
context  of  the  background  medium.  This  derivative  function  can  be  calculated  by 
using  the  ratio  of  two  migration  outputs.  This  ratio  yields  the  stationary  value  of  the 
expression  for  the  derivative  function— exactly  the  result  we  need.  This  is  the  same 
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method  as  is  used  in  Bleistein,  et  al.,  (1987)  to  compute  the  cosine  of  the  incident 
angle  at  a  reflector. 

The  efficiency  of  a  velocity  analysis  approach  surely  depends  on  the  choice  of  mi¬ 
gration  algorithm  used.  In  general,  integral-type  migration  approaches,  such  as  Kirch- 
hoff  or  Gaussian  beam,  are  preferable,  because  those  methods  can  be  implemented 
either  in  common-shot  gathers  or  in  common-offset  gathers,  and  have  the  flexibility  to 
image  the  targeted  structures  in  which  velocity  is  being  estimated.  Besides  efficiency, 
accuracy  should  be  considered  as  well.  For  imaging  complex  structures,  a  migration 
algorithm  should  be  designed  to  handle  turning  waves  and  caustic  regions. 

In  the  Kirchhoff  integral,  traveltime  calculation  by  ray  tracing  or  finite  differencing 
dominates  the  total  cost.  Finite  differencing  calculates  only  the  first-arrival  travel¬ 
time,  so  this  approach  works  fast  but  fails  in  calculation  of  major-energy  traveltimes 
in  caustic  regions  (Geoltrain  and  Brae,  1993).  The  paraxial  ray  method,  in  contrast, 
shows  advantages  in  handling  multivalued  traveltimes  and  caustics.  This  method  uses 
information  from  the  standard  dynamic  ray-tracing  method  to  extrapolate  traveltimes 
and  ray  amplitudes  at  receivers  in  the  vicinity  of  a  central  ray  (Beydoun  and  Keho, 
1987).  However,  the  traveltime  calculation  by  the  paraxial  ray  tracer  is  generally 
more  costly  than  that  by  finite  differencing.  Here,  finite  differencing  is  initially  used 
in  the  migration  implementation  for  velocity  analysis  in  the  area  of  the  Marmousi 
model  where  the  first  arrivals  carry  the  major  energy.  In  the  central  bottom  parts, 
a  paraxial  ray  tracing  algorithm  is  used  to  implement  Kirchhoff  migration.  In  this 
approach,  the  traveltime  corresponding  to  the  major  energy  is  chosen  when  multiple 
arrivals  exist. 


METHODOLOGY 

If  an  incorrect  velocity  is  used  in  prestack  migration,  the  imaged  depths  from 
different  offsets  in  a  common-image-gather(CIG)  will  differ  from  each  other.  A  CIG 
is  a  gather  in  which  migrated  traces  have  the  same  lateral  image  location.  In  this 
situation,  a  residual  moveout  is  observed  in  migrated  data.  The  principle  of  velocity 
analysis  is  to  correct  the  velocity  so  that  the  imaged  depths  at  each  common-image- 
gather  are  close  to  each  other.  For  this  purpose,  one  needs  a  quantitative  relationship 
between  the  residual  moveout  and  the  velocity  error.  Here,  we  present  a  representa¬ 
tion  for  residual  moveout  derived  by  using  a  perturbation  method.  This  representation 
is  valid  for  any  offset,  reflector  dip  and  velocity  distribution. 

Perturbation  formula  for  velocity  estimation 

We  denote  by  a;  a  2-D  vector,  x  =  (x,  z ).  Let  xa  be  the  source  position  and  xr  be 
the  receiver  position  on  the  surface.  For  any  point  X  below  the  surface,  ra(xa,x )  or 
rr(x,xr),  respectively,  denote  traveltimes  from  xs  to  x,  or  x  to  xr. 

Suppose  we  know  the  total  reflection  travetime  function  t(y,h)  (therefore,  dt/dy) 
that  depends  on  midpoint  y  and  half-offset  h.  Given  a  velocity  function  v(x,  2),  then, 
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for  each  h,  the  reflector  is  determined  by 


x)  +  Tr(x,xr)  =  t(y,h), 

(1) 

dr,  drr  dt 
dy  +  dy  dy' 

(2) 

In  a  common  image  gather,  the  imaged  depth  z  can  be  determined  as  a  function  of 
h.  If  the  migration  velocity  equals  the  true  velocity,  then  the  imaged  depth  z  will  be 
independent  of  offset  h\  otherwise — for  incorrect  velocity — 2  varies  with  offset  h.  Con¬ 
sequently,  the  imaged  depths  in  common  image  gathers  (CIGs)  provide  information 
on  velocity  distribution. 

Equations  (1)  and  (2)  are  valid  even  if  the  migration  velocity  is  wrong.  The 
dependency  of  traveltime  on  velocity  in  equations  (1)  and  (2)  implies  that  this  equa¬ 
tion  system  displays  a  general  relationship  between  the  imaged  depth  and  migration 
velocity.  However,  this  equation  system  is  nonlinear,  making  it  difficult  to  solve  di¬ 
rectly  for  velocity.  Here,  we  use  a  mathematical  tool— perturbation — to  linearize  this 
equation  system  by  considering  for  all  perturbations  in  model  parameters. 

Suppose  that  the  velocity  distribution  v  is  characterized  by  a  parameter  or  a  family 
of  parameters,  A, 

v  =  v{x\  A). 

For  example,  when  v(x ;  A)  =  v0  +  ax  +  bz,  A  is  any  set  of  one  to  three  parameters 
chosen  from  vq,  a,  and  b.  Thus,  the  problem  of  velocity  estimation  becomes  a  problem 
of  parameter  estimation.  If  there  is  a  small  perturbation  <5 A  =  A*  —  A  between  the 
true  parameter  and  the  parameter  used  in  migration,  then  the  imaged  depth  will  have 
a  corresponding  perturbation 

6z  =  g(x,h)8\.  (3) 

The  derivative  function  g  can  be  determined  based  on  equations  (1)  and  (2).  To 
simplify  the  derivation,  we  suppose  that  A  is  just  a  single  parameter  at  first. 

For  a  fixed  image  location  x,  we  differentiate  equation  (1)  with  respect  to  A.  Note 
that  y  and  z  are  functions  of  A;  then 


drs  8tt 
dy  dy 


dr,  drT 
d\  +  d\ 


dr s  drr 
dz  +  dz 


dz 

dX 


dt  dy 
dydX 


(4) 


By  using  equation  (2),  the  first  term  of  the  left  side  in  equation  (4)  is  balanced  by 
the  right-hand  term.  Therefore, 


drs  drT  dz  _  dr3  drr 
dz  dz  dX  dX  dX ' 


(5) 


Let  $s  or  9r  be  the  angle  between  the  raypath  from  the  source  or  the  receiver,  and 
the  vertical  at  x ;  then 

drs  _  cosflj  8tt  _  cos0r 

dz  v(x]X)'  dz  v^iX)’ 
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By  using  the  above  equation,  equation  (5)  is  rewritten  as 


cos  6S  +  cos  6a  dz 


drs 
'  dX 


drr 

~dX' 


v(x;X)  dX 

Thus  the  derivative  of  imaged  depth  with  respect  to  A  is  found  by 

v(x ;  A) 


g(x,  h)  =  - 


dr,  drT 
dX  dX 


cos  03  +  cos  6r 


(7) 

(8) 


The  function  g  characterizes  the  relationship  between  the  imaged  depth  and  the 
migration  velocity  in  a  general  medium  context.  The  computation  of  this  function 
will  result  in  a  new  migration  velocity  analysis  method,  compared  to  conventional 
ones  based  on  hyperbolic  residual  moveout. 


The  velocity  estimation  based  on  equation  (3)  has  no  limitation  in  velocity  dis¬ 
tribution  and  reflection  geometry  as  long  as  the  perturbation  is  small.  When  the 
velocity  distribution  is  characterized  by  multiple  parameters,  A  =  (Ai,  A2,  ...,  An)r, 
the  imaged-depth  perturbation  will  depend  on  the  perturbations  of  all  these  param¬ 
eters.  Therefore,  equation  (3)  will  be  modified  to 


n  dz  " 

Sz(x,h)  =  oT-<5A,-  =  J29i(x,h)6Xi, 

1=1  OAi  i= 1 

where  each  derivative  function  is  calculated  by 


gi(x,h)  =  - 


dr,  drr 
dXi  +  dXi 


v(x;  A) 

cos  63  +  cos  9r ' 


(9) 


(10) 


Calculation  of  the  function  g 

The  function  g(y,h),  which  involves  the  derivatives  of  traveltimes  with  respect  to 
the  parameter  A,  can  be  calculated  by 


where  L  is  the  raypath  from  the  source  (or  receiver)  to  the  image  point  x. 

For  each  source  or  receiver,  drjdX  can  be  determined  from  equation  (11).  There¬ 
fore,  given  an  image  point  x  and  a  specular  source-receiver  pair  xs  and  xr,  one  can 
calculate  g  from  formula  (8).  However,  there  is  no  explicit  formula  to  represent  the 
specular  source-receiver  pair  from  the  image  point  for  a  complex  medium.  To  solve 
this  problem,  we  use  the  Kirchhoff  integral  to  calculate  g.  In  the  KirchhofF summation, 
we  calculate  two  migration  outputs  which  have  the  same  phase  but  different  ampli¬ 
tudes.  One  uses  the  original  amplitude;  the  other  one  uses  the  original  amplitude 
multiplied  by  the  quantity  g.  Thus,  the  ratio  of  the  amplitudes  of  these  two  outputs 
will  evaluate  g  at  the  specular  source-receiver  position  without  requiring  knowledge  of 
the  specular  source-receiver  pair.  This  technique,  based  on  stationary-phase  analysis, 
is  the  same  as  was  used  to  determine  the  angle  of  reflection  in  Kirchhoff  inversion 
(Bleistein  et  al.,  1987). 
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Characterization  of  velocity  distribution 

Although  equation  (9)  holds  for  any  velocity  distribution,  the  solution  will  be 
underdetermined  and  unstable  if  too  many  unknown  parameters  are  involved.  Con¬ 
sequently,  it  is  essential  to  characterize  the  velocity  distribution  by  choosing  only 
a  few  appropriate  parameters.  Conventionally,  one  assumes  that  a  velocity  model 
consists  of  the  construction  of  the  macro-model  (constant  velocities  and  velocity  in¬ 
terfaces).  The  interfaces  divide  the  whole  model  into  a  number  of  blocks  or  layers 
(shown  in  Figure  1). 


Here,  we  replace  constant  velocity  in  each  block  by  a  linear  function  that  is  char¬ 
acterized  by  three  parameters: 


Ai  +  A2  (z  —  Zq)  +  \z(x  —  £0), 

where  ( xq,zo )  is  a  reference  point.  Thus,  the  velocity  distribution  is  written  in  the 
form 

v(x,  z)  =  «o(ar,  z)  +  Ax  +  \2(z  -  zQ)  +  A3(a;  -  x0),  (12) 

where  vq  is  a  background  velocity.  An  imaged  depth  depends  on  only  the  velocity 
above  it,  except  for  turning  rays.  Therefore,  use  of  a  recursive  algorithm  (layer 
stripping)  is  possible  to  determine  velocity  in  an  individual  block.  We  start  from 
the  block  nearest  surface.  For  example,  we  choose  the  top  left  block  in  Figure  1. 
In  each  block,  iteration  is  used  to  calculate  velocity  parameters.  Given  an  initial 
guess  for  A,-’s,  prestack  depth  migration  is  implemented  to  obtain  imaged  depths  and 
gi(x,h )  in  equation  (10)  for  common  image  gathers  that  will  give  a  correction  of  the 
parameters.  Then  by  using  the  updated  parameters  as  an  initial  guess,  we  correct 
the  velocity  again  until  convergence  is  achieved.  After  velocity  analysis  in  one  block, 
we  migrate  data  with  the  corrected  velocity,  and  pick  the  velocity  interface  from  the 


5 


Liu  and  Bleistein 


Velocity  analysis 


imaged  structures.  When  velocity  and  velocity  interface  determining  is  finished  in  one 
block,  we  will  repeat  the  same  procedure  to  the  next  block  that  is  located  directly 
below  the  finished  block. 

The  layer-stripping  procedure  for  velocity  analysis  can  be  stated  as  follows: 

•  begin  with  the  first  block 

1 .  estimate  velocity  parameters  iteratively 

(a)  migrate  with  an  initial  guess  of  velocity; 

(b)  sort  the  migrated  data  into  common  image  gathers; 

(c)  measure  imaged  depths  and  evaluate  the  derivative  function; 

(d)  update  velocity  by  using  the  perturbation  formula; 

2.  image  velocity  interface  by  using  corrected  velocity 

•  repeat  step  1  and  2  for  next  block  (recursion) 

When  velocity  and  velocity  gradients  are  estimated  simultaneously  in  a  given 
block,  the  iteration  tends  to  be  unstable.  To  overcome  this  difficulty,  it  is  preferable 
to  estimate  velocity  first.  This  estimation  will  yield  an  average  velocity  and  give  a 
better  initial  guess  for  the  velocity  gradients  in  this  block. 

COMPUTER  IMPLEMENTATION 

The  Marmousi  data  set  is  generated  by  using  a  two-dimensional  acoustic  finite- 
difference  modeling  program.  The  model  contains  many  reflectors,  steep  dips,  and 
strong  velocity  variations  in  both  lateral  and  vertical  directions  (with  a  minimum 
velocity  of  1500  m/s  and  a  maximum  velocity  of  5500  m/s),  shown  in  Figure  2.  The 
data  set  consists  of  240  shots  with  96  traces  per  shot.  The  initial  offset  is  200  m;  both 
the  shot  and  receiver  spacings  are  25  m.  The  first  shot  is  at  lateral  position  3000 
m.  Nineteen  common-offset  data  gathers  with  offsets  ranging  from  200  m  to  2000  m 
are  used  for  velocity  analysis.  The  selected  offsets  range  from  200  m  to  2000  m  with 
spacing  100  m.  The  minimum-offset  gather  is  shown  in  Figure  3. 

During  the  velocity  analysis  process,  we  assume  that  the  velocity  field  is  a  macro 
model  and  that  the  velocity  distribution  is  a  linear  function  in  each  block.  Veloc¬ 
ity  analysis  results  surely  depend  the  choice  of  migration  algorithm  used.  There  are 
two  commonly  used  approaches  to  calculate  traveltimes  in  Kirchhoff  migration:  finite 
differencing  and  ray  tracing.  Compared  to  ray  tracing,  the  finite  differencing  ap¬ 
proach  is  easier  to  code  and  more  efficient  to  implement,  but  fails  to  correctly  image 
complicated  structures  when  multiple  arrivals  exist  (Geoltrain  and  Brae,  1993). 

A  finite  differencing  traveltime  solver  is  initially  used  in  the  migration  implemen¬ 
tation  for  velocity  analysis  in  the  area  where  the  first  arrivals  carry  the  major  energy. 
The  estimated  velocity  model  is  shown  in  Figure  4.  In  the  central  bottom  parts, 
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a  paraxial  ray  tracing  algorithm  is  used  to  implement  Kirchhoff  migration.  In  this 
approach,  the  traveltime  corresponding  to  the  major  energy  is  chosen  when  multiple 
arrivals  exist. 

Using  the  Kirchhoff  migration  algorithm  by  paraxial  ray  tracing,  velocity  analysis 
is  done  through  the  central  bottom  parts  of  the  Marmousi  model.  The  updated 
velocity  model,  shown  in  Figure  5,  consists  of  19  blocks.  In  each  block,  the  velocity 
distribution  is  a  constant  or  linear  function  of  the  depth.  Comparison  of  the  estimated 
velocity  model  with  the  true  one  at  three  lateral  locations  is  shown  Figure  6.  One  can 
see  that  the  estimated  velocity  matches  the  true  model  well  except  for  thin  layers. 
In  fact,  from  sensitivity  analysis  (Liu,  1995),  velocities  in  these  thin  layers  cannot  be 
determined  well.  The  stacked  migration  section  using  the  velocity  model  in  Figure  5 
is  shown  in  Figure  7.  Compared  to  the  migration  result  using  the  true  velocity 
model,  shown  in  Figure  8,  Figure  7  gives  an  acceptable  structural  image  even  in  the 
central  bottom  parts,  which  indicates  the  capability  of  this  migration  velocity  analysis 
approach  for  handling  complex  structures.  The  subsurfaces  are  well  imaged  except 
for  some  detailed  features  in  the  central  bottom  parts.  Some  blurry  image  in  the 
central  bottom  parts  may  be  caused  by  missing  high-velocity  zones.  Figure  2  shows 
that  the  true  velocity  model  contains  several  small  high-velocity  zones  in  the  central 
parts.  These  high-velocity  zones  do  not  appear  in  Figure  5  because  of  the  limitations 
of  velocity  analysis  and  the  resolution  of  migration  imaging. 

Selected  common  image  gathers  from  migrated  data  using  the  estimated  model 
are  shown  in  Figures  9,  10  and  11,  representing  the  left,  central  and  right  parts  of 
the  model  respectively.  The  alignment  of  reflections  in  Figures  9,  11  and  the  upper 
part  of  Figure  10  indicates  the  correctness  of  the  estimated  velocity  in  these  areas. 
The  bottom  part  of  Figure  10  shows  incoherent  signals  that  limit  the  accuracy  of 
velocity  estimation  in  this  particular  area.  Notice  that,  at  the  same  common-image 
gathers  (see  Figure  12),  migrated  data  using  the  true  velocity  also  contain  obvious 
incoherency  in  residual  moveouts,  although  the  stacked  section  shows  a  good  image. 
This  example  demonstrates  that  for  an  extremely  complex  structure  it  is  difficult  to 
identify  the  correct  velocity  model  based  on  the  criterion  of  kinematic  coherence. 

CONCLUSION 

Imaging  complex  structures  (such  as  the  Marmousi  data)  requires  effective  prestack 
depth  migration  algorithms  as  well  as  advanced  velocity  analysis  techniques.  The  per¬ 
turbation  method  in  this  paper  provides  a  useful  tool  for  updating  a  velocity  model  by 
matching  a  criterion  based  on  kinematic  coherence  of  the  prestack  migrated  images, 
which  is  one  of  the  key  elements  in  velocity  model  determination  (Versteeg,  1994).  In 
order  to  obtain  a  more  satisfactory  velocity  analysis  result  for  a  complicated  model, 
as  Versteeg  concluded,  one  also  needs  related  geologic  information  as  constraints. 
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Fig.  3.  The  minimum-offset  Marmousi  data.  The  offset  is  200  meters. 
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Fig.  4.  The  estimated  velocity  model  using  the  first-arrival  operater. 
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Fig.  5.  The  updated  velocity  model  using  the  paraxial  ray  tracer. 
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Fig.  6.  Comparison  of  the  true  velocity  model  with  the  estimated  one.  The  dark 
curve  denotes  the  true  velocity  and  the  gray  denotes  the  estimated  velocity.  The  top 
figure  is  at  location  x=8  km;  the  middle,  at  location  x=6  km;  the  bottom,  at  location 
x=4  km. 
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Fig.  7.  19-offset  stacked  migration  output  for  the  Marmousi  data.  The  input 

velocity  is  one  in  Figure  6. 


4000 


Midpoint  (m) 

5000  6Q,00  70,0Q 


8000 


^  1000 


Q. 

Q  2000 


3000 


Fig.  8.  19-offset  stacked  migration  output  for  the  Marmousi  data.  The  true  velocity 

is  used. 
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Gathers 


Ten  common  image  gathers.  19  offsets  in  each  CIG.  The  image  location 
ranges  from  4  km  to  4.25  km. 


Fig.  10.  Ten  common  image  gathers.  19  offsets  in  each  CIG.  The  image  location 

ranges  from  6  km  to  6.25  km. 


13 


Fig.  12.  Ten  common  image  gathers  from  migrated  data  using  the  true  velocity. 
The  image  location  ranges  from  6  km  to  6.25  km. 
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